replication_txt
Isabella Muratore
November 26, 2017
library(ggplot2)
library(dplyr)## Warning: package 'dplyr' was built under R version 3.4.2##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(plyr)## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizelibrary(vegan)## Warning: package 'vegan' was built under R version 3.4.2## Loading required package: permute## Warning: package 'permute' was built under R version 3.4.2## Loading required package: lattice## This is vegan 2.4-4library(reshape2)
library(VennDiagram)## Warning: package 'VennDiagram' was built under R version 3.4.2## Loading required package: grid## Loading required package: futile.logger## Warning: package 'futile.logger' was built under R version 3.4.2library(knitr)## Warning: package 'knitr' was built under R version 3.4.3Replication of figures 1, 3, and 4 (and Table 1) of “Gene expression patterns associated with caste and reproductive status in ants: worker-specific genes are more derived than queen-specific ones” by Feldmeyer, Elsner, and Foitzik, 2013
This paper involved sequencing total RNA from four different phenotypes of female ants that occur in colonies of Temnothorax longispinosus
Figure 1a, NMDS plot
as published
In this NMDS plot, all the genes that are significantly differentially expressed (p < 0.05 and fold change > 2) are plotted in four dimensions based on the four expression values each gene has for the four categories sampled, in a manner that should minimize “stress”, or how hard it is for the compressed 2D plot to accurately depict the distances between the points. However, this NMDS particular plot does not result in finding a real minimum. The algorithm finds a somewhat low stress value to settle at but does not establish it as a minimum.
The purpose of this figure is to illustrate any distinctive clusters of genes that appear, indicating differences in gene expression patterns associated with each phenotype. The distance between the red labels, which are actually calculated as “species scores”, or weighted averages for all the expression values belonging to that category, illustrate roughly how different the sampled females are from each other under this mode of analysis.
#filtering for genes that have p<0.05 and fc>2 for at least one of the six different pairwise comparisons. These genes are referred to as "DEGs", differentially expressed genes
b <- filter(a, queenfertile_p<0.05 & abs(queenfertile_fold)>2 | queeninfertile_p<0.05 & abs(queeninfertile_fold)>2 | queenforager_p<0.05 & abs(queenforager_fold)>2 | fertileinfertile_p<0.05 & abs(fertileinfertile_fold)>2 | fertileforager_p<0.05 & abs(fertileforager_fold)>2 | infertileforager_p<0.05 & abs(infertileforager_fold)>2)
#selecting just the expression values
s <- mutate(b) %>% select(queen, fertile, infertile, forager)
#performing nmds
#s.mds <- metaMDS(s)
#s.mds
#plotting using points
#plot(s.mds, type = "p")
#My NMDS plot is almost exactly the same as the published one, but mine is rotated in the opposite orientation
#plotting using text to show that the red labels match up to their positions in the published figure.
#plot(s.mds, type = "t")When I filtered for genes that had p<0.05 and fc>2 for at least one of the six pairwise comparisons, I obtained an NMDS plot almost identical to the published version, except that it was rotated by a multiple of 90 degrees from the published orientation. This does not reflect any meaningful difference.
Figure 1b, Venn diagram
as published
This Venn diagram is supposed to show the subsets of diffentially expressed genes that were expressed in the four different categories and how these lists overlap.
The vegan package I (and the authors) used for the diagram outputs to a file, not directly to R, so I am embedding the resulting image.
# SECOND CLOSEST
q_nonzero <- filter(a, queen != 0 & queenfertile_p<0.05 & abs(queenfertile_fold)>2 & queeninfertile_p<0.05 & abs(queeninfertile_fold)>2 & queenforager_p<0.05 & abs(queenforager_fold)>2)
q_nonzero <- mutate(q_nonzero) %>% select(ID)
fert_nonzero <- filter(a, fertile != 0 & queenfertile_p<0.05 & abs(queenfertile_fold)>2 & fertileinfertile_p<0.05 & abs(fertileinfertile_fold)>2 & fertileforager_p<0.05 & abs(fertileforager_fold)>2)
fert_nonzero <- mutate(fert_nonzero) %>% select(ID)
infert_nonzero <- filter(a, infertile != 0 & fertileinfertile_p<0.05 & abs(fertileinfertile_fold)>2 & queeninfertile_p<0.05 & abs(queeninfertile_fold)>2 & infertileforager_p<0.05 & abs(infertileforager_fold)>2)
infert_nonzero <- mutate(infert_nonzero) %>% select(ID)
for_nonzero <- filter(a, forager != 0 & infertileforager_p<0.05 & abs(infertileforager_fold)>2 & fertileforager_p<0.05 & abs(fertileforager_fold)>2 & queenforager_p<0.05 & abs(queenforager_fold)>2)
for_nonzero <- mutate(for_nonzero) %>% select(ID)
####
#CLOSEST
q_nonzero <- filter(a, queen != 0 & queenfertile_p<0.05 & queeninfertile_p<0.05 & queenforager_p<0.05)
q_nonzero <- mutate(q_nonzero) %>% select(ID)
fert_nonzero <- filter(a, fertile != 0 & queenfertile_p<0.05 & fertileinfertile_p<0.05 & fertileforager_p<0.05)
fert_nonzero <- mutate(fert_nonzero) %>% select(ID)
infert_nonzero <- filter(a, infertile != 0 & fertileinfertile_p<0.05 & queeninfertile_p<0.05 & infertileforager_p<0.05)
infert_nonzero <- mutate(infert_nonzero) %>% select(ID)
for_nonzero <- filter(a, forager != 0 & infertileforager_p<0.05 & fertileforager_p<0.05 & queenforager_p<0.05)
for_nonzero <- mutate(for_nonzero) %>% select(ID)
venn.diagram(list(Queen=q_nonzero$ID, Forager=for_nonzero$ID, Fertile=fert_nonzero$ID, Infertile=infert_nonzero$ID), "replication_venn.tiff", fill = c("blue", "red", "yellow", "green"), cex=3,cat.cex=1.8) ## [1] 1
Figure 3a, pairwise pie
as published
This pie should show the subsets of differentially expressed genes that are statistically different based on p value and fold change between the different pairwise comparisons that can be made.
qfert <- filter(a, queenfertile_p<0.05 & abs(queenfertile_fold)>2)
qfert <- mutate(qfert) %>% select(ID)
qinfert <- filter(a, queeninfertile_p<0.05 & abs(queeninfertile_fold)>2)
qinfert <- mutate(qinfert) %>% select(ID)
qfor <- filter(a, queenforager_p<0.05 & abs(queenforager_fold)>2)
qfor <- mutate(qfor) %>% select(ID)
fertinfert <- filter(a, fertileinfertile_p<0.05 & abs(fertileinfertile_fold)>2)
fertinfert <- mutate(fertinfert) %>% select(ID)
fertfor <- filter(a, fertileforager_p<0.05 & abs(fertileforager_fold)>2)
fertfor <- mutate(fertfor) %>% select(ID)
infertfor <- filter(a, infertileforager_p<0.05 & abs(infertileforager_fold)>2)
infertfor <- mutate(infertfor) %>% select(ID)
slices <- c(nrow(qfert), nrow(qinfert), nrow(qfor), nrow(fertinfert), nrow(fertfor), nrow(infertfor))
lbls <- c("q vs. fer", "q vs. infer", "q vs. for", "fer vs. inf", "fer vs. for", "inf vs. for")
lbls <- paste(lbls, slices) # add percents to labels
pie(slices,labels = lbls, col=rainbow(length(lbls)),
main="Pie Chart of Pairwise Differences")
Figure 3b, up-regulated pie
as published
This pie should show the subsets of differentially expressed genes that are uniquely upregulated in one category relative to all the others. These would suggest that these genes are fundementally important for creating the unique phenotype in question, probably being necessarily for developping certain structures needed for its behavior and role.
#b <- filter(a, queenfertile_p<0.05 & abs(queenfertile_fold)>2 | queeninfertile_p<0.05 & abs(queeninfertile_fold)>2 | queenforager_p<0.05 & abs(queenforager_fold)>2 | fertileinfertile_p<0.05 & abs(fertileinfertile_fold)>2 | fertileforager_p<0.05 & abs(fertileforager_fold)>2 | infertileforager_p<0.05 & abs(infertileforager_fold)>2)
#q_up <- filter(b, (queen > forager) & (queen > fertile) & (queen > infertile))
#fert_up <- filter(b, (fertile > forager) & (fertile > queen) & (fertile > infertile))
#infert_up <- filter(b, (infertile > forager) & (infertile > queen) & (infertile > fertile))
#for_up <- filter(b, (forager > fertile) & (forager > queen) & (forager > infertile))
###
#q_up <- filter(a, queenfertile_p<0.05 & queenfertile_fold>0 & queeninfertile_p<0.05 & queeninfertile_fold>0 & queenforager_p<0.05 & queenforager_fold>0)
#fert_up <- filter(a, queenfertile_p<0.05 & queenfertile_fold < 0 & fertileinfertile_p<0.05 & fertileinfertile_fold>0 & fertileforager_p<0.05 & fertileforager_fold>0)
#infert_up <- filter(a, fertileinfertile_p<0.05 & fertileinfertile_fold < 0 & queeninfertile_p<0.05 & queeninfertile_fold < 0 & infertileforager_p<0.05 & infertileforager_fold>0)
#for_up <- filter(a, infertileforager_p<0.05 & infertileforager_fold < 0 & fertileforager_p<0.05 & fertileforager_fold < 0 & queenforager_p<0.05 & queenforager_fold < 0)
###
q_up <- filter(a, queenfertile_p<0.05 & queenfertile_fold>2 & queeninfertile_p<0.05 & queeninfertile_fold>2 & queenforager_p<0.05 & queenforager_fold>2)
fert_up <- filter(a, queenfertile_p<0.05 & queenfertile_fold < -2 & fertileinfertile_p<0.05 & fertileinfertile_fold>2 & fertileforager_p<0.05 & fertileforager_fold>2)
infert_up <- filter(a, fertileinfertile_p<0.05 & fertileinfertile_fold < -2 & queeninfertile_p<0.05 & queeninfertile_fold < -2 & infertileforager_p<0.05 & infertileforager_fold>2)
for_up <- filter(a, infertileforager_p<0.05 & infertileforager_fold < -2 & fertileforager_p<0.05 & fertileforager_fold < -2 & queenforager_p<0.05 & queenforager_fold < -2)
slices <- c(nrow(q_up), nrow(infert_up), nrow(for_up), nrow(fert_up))
lbls <- c("q", "inf", "for", "fer")
lbls <- paste(lbls, slices) # add percents to labels
pie(slices,labels = lbls, col=rainbow(length(lbls)),
main="Pie Chart of Up-Regulated Genes")
Table 1, up-regulation table
as published
This table shows “Number of significantly up-regulated genes (FDR-P < 0.05; fold change >2) in the focal caste in comparison withthe other castes and the number of pair-specific genes (genes up-regulated in only this specific comparison) with the corresponding percentage in parentheses.”
Unlike the previous depictions of numbers of differentially expressed genes between phenotypes, this table illustrates both “directions” of the comparison to differentiate between genes upregulated in one of the two compared phenotypes vs. those upregulated in the other.
# I am filtering for p<0.05 and fc>2 for each pairwise comparison in both "directions", then editing the set down to only contain the column list of gene IDs.
qfert <- filter(a, queenfertile_p<0.05 & queenfertile_fold>2)
qfert <- mutate(qfert) %>% select(ID)
fertq <- filter(a, queenfertile_p<0.05 & queenfertile_fold < -2)
fertq<- mutate(fertq) %>% select(ID)
qinfert <- filter(a, queeninfertile_p<0.05 & queeninfertile_fold>2)
qinfert <- mutate(qinfert) %>% select(ID)
infertq <- filter(a, queeninfertile_p<0.05 & queeninfertile_fold < -2)
infertq <- mutate(infertq) %>% select(ID)
qfor <- filter(a, queenforager_p<0.05 & queenforager_fold>2)
qfor <- mutate(qfor) %>% select(ID)
forq <- filter(a, queenforager_p<0.05 & queenforager_fold < -2)
forq <- mutate(forq) %>% select(ID)
fertinfert <- filter(a, fertileinfertile_p<0.05 & fertileinfertile_fold>2)
fertinfert <- mutate(fertinfert) %>% select(ID)
infertfert <- filter(a, fertileinfertile_p<0.05 & fertileinfertile_fold < -2)
infertfert <- mutate(infertfert) %>% select(ID)
fertfor <- filter(a, fertileforager_p<0.05 & fertileforager_fold>2)
fertfor <- mutate(fertfor) %>% select(ID)
forfert <- filter(a, fertileforager_p<0.05 & fertileforager_fold < -2)
forfert <- mutate(forfert) %>% select(ID)
infertfor <- filter(a, infertileforager_p<0.05 & infertileforager_fold>2)
infertfor <- mutate(infertfor) %>% select(ID)
forinfert <- filter(a, infertileforager_p<0.05 & infertileforager_fold < -2)
forinfert <- mutate(forinfert) %>% select(ID)
#Now I am filtering for the subsets of these gene lists that only appear within that selected comparison and not in any of the others, using dplyr's setdiff function.
q_fert <- dplyr::setdiff(qfert, union(fertq, qinfert, infertq, qfor, forq, fertinfert, infertfert, fertfor, forfert, infertfor, forinfert))
fert_q <- dplyr::setdiff(fertq, union(qfert, qinfert, infertq, qfor, forq, fertinfert, infertfert, fertfor, forfert, infertfor, forinfert))
q_infert <- dplyr::setdiff(qinfert, union(qfert, fertq, infertq, qfor, forq, fertinfert, infertfert, fertfor, forfert, infertfor, forinfert))
infert_q <- dplyr::setdiff(infertq, union(qfert, fertq, qinfert, qfor, forq, fertinfert, infertfert, fertfor, forfert, infertfor, forinfert))
q_for <- dplyr::setdiff(qfor, union(qfert, fertq, qinfert, infertq, forq, fertinfert, infertfert, fertfor, forfert, infertfor, forinfert))
for_q <- dplyr::setdiff(forq, union(qfert, fertq, qinfert, infertq, qfor, fertinfert, infertfert, fertfor, forfert, infertfor, forinfert))
fert_infert <- dplyr::setdiff(fertinfert, union(qfert, fertq, qinfert, infertq, qfor, forq, infertfert, fertfor, forfert, infertfor, forinfert))
infert_fert <- dplyr::setdiff(infertfert, union(qfert, fertq, qinfert, infertq, qfor, forq, fertinfert, fertfor, forfert, infertfor, forinfert))
fert_for <- dplyr::setdiff(fertfor, union(qfert, fertq, qinfert, infertq, qfor, forq, fertinfert, infertfert, forfert, infertfor, forinfert))
for_fert <- dplyr::setdiff(forfert, union(qfert, fertq, qinfert, infertq, qfor, forq, fertinfert, infertfert, fertfor, infertfor, forinfert))
infert_for <- dplyr::setdiff(infertfor, union(qfert, fertq, qinfert, infertq, qfor, forq, fertinfert, infertfert, fertfor, forfert, forinfert))
for_infert <- dplyr::setdiff(forinfert, union(qfert, fertq, qinfert, infertq, qfor, forq, fertinfert, infertfert, fertfor, forfert, infertfor))
#This will be the row in the table showing the total genes that filtered as signficant for each pairwise comparison
#Simultaneously, I am reording the gene sets to match the order in the published table
Total <- c(nrow(qfert), nrow(qinfert), nrow(qfor), nrow(fertq), nrow(fertinfert), nrow(fertfor), nrow(infertfert), nrow(infertq), nrow(infertfor),nrow(forq), nrow(forfert), nrow(forinfert))
#This will be the row in the table showing just the genes uniquely filtering into each pairwise comparison
Pair_specific <- c(nrow(q_fert), nrow(q_infert), nrow(q_for), nrow(fert_q), nrow(fert_infert), nrow(fert_for), nrow(infert_fert), nrow(infert_q), nrow(infert_for),nrow(for_q), nrow(for_fert), nrow(for_infert))
#This will be the the row showing what percent of each gene set is unique to that comparison by dividing one vector by the other
Pair_specific_percent <- Pair_specific/Total
#These vectors match the labels used in the published table
Focal_caste <- c("Queen", "", "", "Fertile worker", "", "", "Infertile worker", "", "", "Forager", "", "")
Comparison <- c("fer", "inf", "fer", "q", "inf", "for", "q", "fer", "for", "q", "fer", "inf")
table1 <- data.frame(Focal_caste, Comparison, Total, Pair_specific, Pair_specific_percent)
#kable function produces a pretty table
knitr::kable(table1)| Focal_caste | Comparison | Total | Pair_specific | Pair_specific_percent |
|---|---|---|---|---|
| Queen | fer | 1055 | 566 | 0.5364929 |
| inf | 876 | 376 | 0.4292237 | |
| fer | 968 | 462 | 0.4772727 | |
| Fertile worker | q | 1728 | 1717 | 0.9936343 |
| inf | 273 | 193 | 0.7069597 | |
| for | 236 | 170 | 0.7203390 | |
| Infertile worker | q | 306 | 118 | 0.3856209 |
| fer | 2043 | 667 | 0.3264807 | |
| for | 194 | 127 | 0.6546392 | |
| Forager | q | 1753 | 416 | 0.2373075 |
| fer | 221 | 57 | 0.2579186 | |
| inf | 144 | 77 | 0.5347222 |
Figure 4, vitellogenin bar graph
as published
Here the authors focus in on a small group of differentially expressed genes since they relate heavily to female ant phenotype. Vitellogenenin is needed for producing eggs, so it is very relevant to reproductive status. This graph shows expression of different copies of the vitellogenin protein and the vitellogenin receptor.
I melted my existing data frame to fit it into the format required for a ggplot grouped bar graph.
#selecting the vitellogenin genes that the authors focused on
s <- filter(a, ID=="TlongiContigs_rep_c8852" | ID=="TlongiContigs_rep_c14651" | ID=="TlongiContigs_rep_c6184" | ID=="TlongiContigs_rep_c39314" | ID=="TlongiContigs_rep_c43820")
#reordering to fit paper
levels(s$ID) <- c("TlongiContigs_rep_c8852", "TlongiContigs_rep_c14651", "TlongiContigs_rep_c6184", "TlongiContigs_rep_c39314", "TlongiContigs_rep_c43820")
#the published figure contains significance groupings denoted by letters, so we should check the pairwise comparison p values to see whether the same groupings appear
#for each gene I am extracting the p values to compare them
Vg2 <- filter(s, ID=="TlongiContigs_rep_c8852") %>% select(ID, queenfertile_p, queeninfertile_p, queenforager_p, fertileinfertile_p, fertileforager_p, infertileforager_p)
Vg2## ID queenfertile_p queeninfertile_p queenforager_p
## 1 TlongiContigs_rep_c8852 4.33584e-13 3.09115e-12 0
## fertileinfertile_p fertileforager_p infertileforager_p
## 1 2.08868e-06 9.88682e-07 6.00246e-13#every expression value is significantly different from the others
Vg3 <- filter(s, ID=="TlongiContigs_rep_c14651") %>% select(ID, queenfertile_p, queeninfertile_p, queenforager_p, fertileinfertile_p, fertileforager_p, infertileforager_p)
Vg3## ID queenfertile_p queeninfertile_p queenforager_p
## 1 TlongiContigs_rep_c14651 2.97754e-12 2.53228e-12 5.25948e-12
## fertileinfertile_p fertileforager_p infertileforager_p
## 1 0.2228907 3.90965e-07 6.52258e-13#fertile and infertile expression values are not significantly different from each other; all others are significantly different
VgRec <- filter(s, ID=="TlongiContigs_rep_c6184") %>% select(ID, queenfertile_p, queeninfertile_p, queenforager_p, fertileinfertile_p, fertileforager_p, infertileforager_p)
VgRec## ID queenfertile_p queeninfertile_p queenforager_p
## 1 TlongiContigs_rep_c6184 7.27345e-13 0 0
## fertileinfertile_p fertileforager_p infertileforager_p
## 1 1 0.1786372 1#queen expression value is signficantly different from all others; others are not signficantly different from each other
Vg1 <- filter(s, ID=="TlongiContigs_rep_c39314") %>% select(ID, queenfertile_p, queeninfertile_p, queenforager_p, fertileinfertile_p, fertileforager_p, infertileforager_p)
Vg1## ID queenfertile_p queeninfertile_p queenforager_p
## 1 TlongiContigs_rep_c39314 0.02210215 4.30138e-13 2.79332e-12
## fertileinfertile_p fertileforager_p infertileforager_p
## 1 1.08503e-05 0.000242488 1#infertile and forager expression values are not significantly different from each other; all others are signficantly different
Vg6 <- filter(s, ID=="TlongiContigs_rep_c43820") %>% select(ID, queenfertile_p, queeninfertile_p, queenforager_p, fertileinfertile_p, fertileforager_p, infertileforager_p)
Vg6## ID queenfertile_p queeninfertile_p queenforager_p
## 1 TlongiContigs_rep_c43820 0.001709597 0.8957542 1.71581e-08
## fertileinfertile_p fertileforager_p infertileforager_p
## 1 3.0524e-05 1.71499e-13 6.08174e-06#queen and infertile values are not significantly different from each other; all others are significantly different
#selecting just the lists of gene IDs
s <- mutate(s) %>% select(ID, queen, fertile, infertile, forager)
#reformatting
new_melt <- melt(s, id.vars='ID')
#gene expression is usually plotting by either log transforming the values or using a log scale, which is what the authors did. This is because values can vary by many orders of magnitude.
g <- ggplot(new_melt, aes(ID, value)) + geom_bar(aes(fill = variable),
width = 0.4, position = position_dodge(width=0.5), stat="identity") +
scale_x_discrete(limits=c("TlongiContigs_rep_c8852", "TlongiContigs_rep_c14651", "TlongiContigs_rep_c6184", "TlongiContigs_rep_c39314", "TlongiContigs_rep_c43820"), labels=c("Vg2", "Vg3", "VgRec", "Vg1", "Vg6"), name="") + scale_y_log10(name="log RPKM")
#adding in significance groupings based on analysis of p values
g <- g + annotate("text",x=0.81,y=1550,label="a")+
annotate("text",x=0.94,y=88,label="c")+
annotate("text",x=1.06,y=200,label="b")+
annotate("text",x=1.19,y=19,label="d")
g <- g + annotate("text",x=1.81,y=1167,label="a")+
annotate("text",x=1.94,y=75,label="b")+
annotate("text",x=2.06,y=123,label="b")+
annotate("text",x=2.19,y=11,label="c")
g <- g + annotate("text",x=2.81,y=127,label="a")+
annotate("text",x=2.94,y=21,label="b")+
annotate("text",x=3.06,y=13,label="b")+
annotate("text",x=3.19,y=4.4,label="b")
g <- g + annotate("text",x=3.81,y=104,label="c")+
annotate("text",x=3.94,y=177,label="b")+
annotate("text",x=4.06,y=320,label="a")+
annotate("text",x=4.19,y=290,label="a")
g <- g + annotate("text",x=4.81,y=74,label="b")+
annotate("text",x=4.94,y=157,label="a")+
annotate("text",x=5.06,y=70,label="b")+
annotate("text",x=5.19,y=12,label="c")
g 